#ifndef _HELPER_H_
#define _HELPER_H_

#include "Spline.h"
#include "SplineHelper.h"
#include <iostream>
#include <vector>
#include <utility>
#include <cmath>


typedef double (*fpd)(double);

double Bfunc(const double x){
     return 1.0/(1+25*x*x);
}

auto DataGenerator(fpd f, const double start, const double end, const
int num){
     std::vector<double> v1,v2;
     const double h = (end - start)/num;
     double t = start;
     std::vector<Vec<double,2> > res;
     for (int i = 0 ; i <= num ; i++){
       res.push_back({t,f(t)});
       t += h;
     }
     return res;
}

template<int Order>
void fnplts(const std::vector<Spline<1,Order,ppForm> >& vsp, const std::string
            &_file,const int num = 200){
  const int N = vsp.size();
  if (N == 1)
    fnplt(vsp[0],_file,num);
  else if (N > 1){
    const std::vector<double>& spknots = vsp[0].getknots();
    std::ofstream os(_file);
    double h = (spknots.back() - spknots.front())/(1.0*num);
    std::vector<std::vector<double> > y(N);
    double t = spknots.front();
    os << "t = [";
    for (int i = 0; i <= num ; i++){
      for (int j = 0 ; j < N ; j++){
        Vec<double,1> tmpvec = vsp[j](t);
        y[j].push_back(tmpvec[0]);
      }
      os << t << ",";
      t += h;
    }
    os << "];\n" << "y = [";
    for (int j = 0; j < N ; j++){
      os << "[";
      for (int i = 0; i <= num ; i++){
        os << y[j][i] << ",";
      }
      os << "],\n";
    }
    os << "];\n" << "plot(t,y);\n";
    os.close();
  } 
}


void Bplot(const std::vector<Spline<1,4,ppForm> >& vsp, const std::string
            &_file,const int num = 200){
  fnplts(vsp,_file,num);  
  ofstream os(_file,ios::app);
  os << "x=-1:0.025:1;\n";
  os << "hold on;\n";
  os << "plot(x,1./(1+25.*x.^2),'*');legend('N=6','N=11','N=21','N=41','N=81','exact function');\n";
  os.close();
}


double Berror(const Spline<1,4,ppForm>& sp, const int N){
  const double h = 2.0/N;
  double t = -1.0 + h/2;
  double Max = 0;
  for (int i = 0 ; i < N ; i++){
    Max = std::max(Max,fabs(sp(t)[0]-Bfunc(t)));
    t += h;
  }
  return Max;
}

double Cfunc(const double x){
  return 1.0/(1+x*x);
}

auto Cgenerator1(){
  std::vector<double> x;
  std::vector<double> y;
  x.push_back(-5);
  y.push_back(5.0/338);
  for (int i = -5 ; i <= 5 ; i++){
    x.push_back(i);
    y.push_back(Cfunc(i));
  }
  x.push_back(5);
  y.push_back(-5.0/338);
  return std::make_pair(x,y);  
}

auto Cgenerator2(){
  std::vector<double> x;
  std::vector<double> y;
  x.push_back(-5);
  y.push_back(Cfunc(-5));
  for (int i = -5 ; i < 5 ; i++){
    x.push_back(i+0.5);
    y.push_back(Cfunc(i+0.5));
  }
  x.push_back(5);
  y.push_back(Cfunc(5));
  return std::make_pair(x,y);  
}

template <int Ord>
void Cplot(const Spline<1,Ord,BForm>& Bs, const std::string
                    &_file, const int num = 200){
  fnplt(Bs,_file,num);  
  ofstream os(_file,ios::app);
  os << "\n";
  os << "x=-5:0.1:5;\n";
  os << "hold on;\n";
  os << "plot(x,1./(1+x.^2),'*');legend('B-spline','exact function');\n";
  os.close();
}

template <int Ord>
double Derror(const Spline<1,Ord,BForm>& Bs, const double t){
  return fabs(Bs(t) - Cfunc(t));
}

double Efunc1(const double x){
  double tmp = std::sqrt(fabs(x))-std::sqrt(3-x*x);
  return 2.0*tmp/3;
}

double Efunc2(const double x){
  double tmp = std::sqrt(fabs(x))+std::sqrt(3-x*x);
  return 2.0*tmp/3;
}

auto Egenerator(const int num){
  assert(num%2 == 0);
  std::vector<Vec<double,2> > vdata1;
  std::vector<Vec<double,2> > vdata2;
  vdata1.push_back({0,-2.0*std::sqrt(3)/3});
  vdata2.push_back({0,-2.0*std::sqrt(3)/3});
  const int N2 = (num - 4)/3;
  const int N1 = num/2 - 2 -N2;
  const double h1 = std::sqrt(3)/(N1+1);
  const double h2 = std::sqrt(3)/(N2+1);
  double t = 0;
  for (int i = 0; i < N1 ; i++){
    t += h1;
    double tmp = Efunc1(t);
    vdata1.push_back({t,tmp});
    vdata2.push_back({-t,tmp});
  }
  vdata1.push_back({std::sqrt(3),2.0*std::sqrt(std::sqrt(3))/3});
  vdata2.push_back({-std::sqrt(3),2.0*std::sqrt(std::sqrt(3))/3});
  t = std::sqrt(3);
  for (int i = 0; i < N2 ; i++){
    t -= h2;
    double tmp = Efunc2(t);
    vdata1.push_back({t,tmp});
    vdata2.push_back({-t,tmp});
  }
  vdata1.push_back({0,2.0*std::sqrt(3)/3});
  vdata2.push_back({0,2.0*std::sqrt(3)/3});
  return std::make_pair(vdata1,vdata2);
}

void Eplot(const std::vector<std::pair<Spline<2,4,ppForm>,Spline<2,4,ppForm> > >& vpsp,const std::string
           &_file, const int num = 200){
  ofstream os(_file);
  os << "figure;\n";
  os << "subplot(1,3,1);\n";
  os.close();
  fnplt(vpsp[0].first,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "hold on;\n";
  os.close();
  fnplt(vpsp[0].second,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "title('n=10');\n";
  os << "subplot(1,3,2);\n";
  os.close();
  fnplt(vpsp[1].first,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "hold on;\n";
  os.close();
  fnplt(vpsp[1].second,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "title('n=40');\n";
  os << "subplot(1,3,3);\n";
  os.close();
  fnplt(vpsp[2].first,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "hold on;\n";
  os.close();
  fnplt(vpsp[2].second,_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "title('n=160');\n";
  os.close();
}

void Fplot(const std::vector<Spline<1,4,ppForm> >& vsp,const std::string
           &_file, const int num = 200){
  fnplt(vsp[0],_file,num);
  ofstream os(_file,ios::app);
  os << "\n" << "hold on;\n";
  os.close();
  fnplt(vsp[1],_file,num,1);
  os.open(_file,ios::app);
  os << "\n" << "hold on;\n";
  os.close();
  fnplt(vsp[2],_file,num,1);
}

template<int N>
void Gplot(const std::string &_file, const int num = 200){
  const int size = N + 2;
  const double left = -3;
  const double right = 5;
  const int start = 0;
  std::vector<Spline<1,N+1,ppForm> > table(size*size);
  for(int i = 0; i < size ; i++){
    Spline<1,N+1,ppForm> tmpsp;
    std::vector<Polynomial<N+1,Vec<double,1> > > tmppoly(2);
    Polynomial<N+1,Vec<double,1> > basis;
    basis[N] = Vec<double,1>({-left+i});
    basis[N-1] = Vec<double,1>({-1});
    for (int j = 0 ; j < N - 1 ; j++)
      basis = adjust<N+1,2*N+1>(basis*basis);
    tmppoly[0] = basis;
    tmppoly[1] = Polynomial<N+1,Vec<double,1> >();
    tmpsp.setpolys(tmppoly);
    std::vector<double> tmpknot(3);
    tmpknot[0] = left;
    tmpknot[1] = start + i;
    tmpknot[2] = right;
    tmpsp.setknots(tmpknot);
    table[i] = tmpsp;
  }
  for (int i = 1; i < size ; i++){
    for (int j = i ; j < size ; j++){
      Spline<1,N+1,ppForm> tmpsp = table[(i-1)*size+j] -
      table[(i-1)*size+j-1];
      table[i*size+j] = tmpsp;
    }
  }
  ofstream os(_file);
  os << "figure;\n";
  os << "axis off;\n";
  for (int i = 0 ; i < size ; i++){
    for (int j = i ; j < size ; j++){
      os << "subplot(" << size << "," << size << "," << i+1+j*size <<
      ");\n";
      os.close();
      fnplt(table[i*size+j],_file,num,1);
      os.open(_file,ios::app);
      os << "axis([-3,4,-0.5,3]);\n";
    }    
  }
  os.close();
}



#endif // _HELPER_H_
